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LONG-TERM  GOALS 

My  long  term  goal  is  to  develop  sea  ice  dynamics  models  that  describe  behavior  on  kilometer  scales 
and  larger,  and  to  base  these  models  on  the  smaller  scale  physical  processes  known  to  control  leading, 
rafting,  and  ridging.  I  also  aim  to  implement  these  models  in  the  next- generation  Polar  Ice  Prediction 
System  (PIPS). 

OBJECTIVES 

My  contract  has  the  following  tasks:  (1)  complete  development  of  the  constitutive  law  and  test  it  with 
Od  calculations,  (2)  incorporate  the  model  into  a  2d  numerical  code  and  test  its  performance  by 
comparing  behavior  with  observed  behavior,  (3)  extend  my  ambient  noise  model  by  relating  the  noise¬ 
generating  model  to  ice  mechanical  behavior,  (4)  help  transfer  the  new  models  into  fleet  systems  such 
as  PIPS,  and  (5)  publish  the  results.  The  contract  was  expanded  during  1999  to  accelerate 
development  of  the  ice  dynamics  model  for  PIPS  3.0.  My  tasks  include  adaptation  of  an  implicit  finite 
element  method  and  reformulation  of  the  redistribution  function  in  the  oriented  thickness  distribution 
description. 

APPROACH 

We  have  developed  a  new  anisotropic  plasticity  constitutive  law  to  describe  and  forecast  ice  stress, 
deformation,  lead  direction,  and  ice  condition  at  scales  from  a  few  kilometers  to  hundreds  of 
kilometers  [Coon,  etal.,  1992,  1998;  Pritchard,  1998a].  The  essential  difference  between  this 
anisotropic  plasticity  model  and  previous  isotropic  models  is  that  it  can  describe  the  formation  and 
direction  of  each  new  lead  or  ridge  system  and  track  its  thickness  distribution.  I  am  working  with  M. 
Coon  to  develop  this  constitutive  law. 

Ice  conditions  are  described  by  an  oriented  thickness  distribution  g{h,b )  that  describes  the  fraction  of 
ice  that  has  thickness  h  and  orientation  0  [Pritchard,  1998a].  This  description  is  needed  for  Eulerian 
simulations,  where  we  must  advect  complete  fields  through  the  grid. 

I  have  been  reviewing  numerical  methods  used  in  the  finite  element  codes  to  determine  how  to  adapt 
an  implicit  scheme  to  the  new  ice  dynamics  model.  Several  methods  are  attractive.  They  offer  us  the 
opportunity  to  use  time  steps  representative  of  the  physical  processes,  perhaps  on  the  order  of  a  day. 
The  algorithms  are  more  complicated,  and  each  step  is  slower  than  the  existing  explicit  leap  frog 
scheme.  The  new  methods  will  be  tested  in  a  finite  element  code  that  I  developed  earlier  for  ice 
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dynamics  simulations.  I  am  talking  with  Dr.  H.  L.  Schreyer  of  the  University  of  New  Mexico  to  learn 
how  the  material  science  community  has  adapted  these  methods  to  their  problems. 

I  have  been  working  with  R.  Bourke,  J.  Wilson,  and  W.  D.  Hibler  III  to  modify  the  PIPS  2  ice 
dynamics  model  so  that  it  can  serve  as  an  acoustic  source  model.  My  experience  using  energy 
dissipation  as  a  proxy  variable  to  describe  noise  generated  by  ridging  [Pritchard,  1984,  1990,  1993] 
provides  this  team  with  the  needed  background.  These  results  will  provide  guidance  to  development  of 
an  acoustical  model  for  PIPS  3.0. 

WORK  COMPLETED 

The  anisotropic  elastic  plastic  constitutive  law  and  the  oriented  thickness  distribution  have  been 
formulated  [Coon,  et  al.,  1992,  1998;  Pritchard,  1998a].  I  have  developed  a  code  to  integrate  the 
anisotropic  plasticity  constitutive  law  for  given  deformation  histories  [Pritchard,  1998b].  This  Od 
model,  which  is  similar  to  the  column  models  used  in  large-scale  climate  simulations  will  be  used  to 
define  material  properties,  in  preparation  for  full  2d  simulations. 

The  redistribution  function  has  been  restudied,  and  a  new  formulation  developed.  I  have  extended  the 
formulation  presented  by  Pritchard  [1998a]  to  correct  the  behavior  under  all  ice  conditions.  My  earlier 
function  was  not  adequate  when  thin  ice  in  a  lead  was  lacking. 

I  participated  in  a  PIPS  3.0  Development  Team  meeting  at  the  National  Ice  Center.  This  group 
discussed  concepts  for  the  next-generation  ice  dynamics  model,  and  the  anisotropic  plasticity  law  will 
be  used  in  PIPS  3.0. 

A  finite  element  numerical  method  has  been  identified  as  a  candidate  for  integrating  the  next- 
generation  anisotropic  plasticity  model. 

RESULTS 

Redistribution  Function.  At  each  instant  and  location,  the  plasticity  law  determines  the  stress  state  and 
the  plastic  stretching  D;)  (the  importance  of  this  fact  has  not  always  been  recognized).  The  stress  state 

at  plastic  failure  determines  the  orientation  of  ice  that  is  deforming.  There  can  be  no  stretching  along  a 
lead  or  ridge  unless  the  surrounding  thicker  ice  is  failing.  If  the  plastic  stretching  is  composed  of 
contributions  from  more  than  one  flow  rule,  then  ice  may  be  failing  in  multiple  orientations.  This  may 
happen  when  the  stress  is  on  a  comer  of  the  yield  surface,  or  it  is  on  an  isotropic  part  of  the  yield 
surface. 

The  plastic  stretching  is  then  written  as 


D,.  =  Ip, 

j 


is  the  contribution  to  plastic  stretching  from  the  0,- orientation. 


If  the  plastic  stretching  is  opening,  then  open  water  forms  in  that  orientation  (or  orientations).  If  the 
plastic  stretching  is  closing,  thin  ice  in  the  lead  having  that  orientation  will  ridge,  or  if  all  of  this  thin 
ice  has  been  consumed,  then  ice  having  all  orientations  (isotropic  ice)  will  ridge.  This  is  the 
fundamental  extension  to  the  redistribution  function  introduced  by  Pritchard  [1998a].  When  isotropic 
ice  is  being  deformed,  the  behavior  can  follow  the  isotropic  thickness  distribution  theory  introduced  by 
Thorndike,  et  al.  [1975].  For  the  oriented  thickness  distribution  g(h,&) ,  the  isotropic  fractions  are 
defined  by 

i(h)  =  min  g(h,e ) . 

e 

We  assume  that  only  the  thinner  (or  thinnest)  fractions  participate  in  the  ridging.  A  manuscript 
describing  the  new  redistribution  function  is  now  being  prepared. 

For  each  deforming  orientation,  the  plastic  stretching  should  be  expressed  in  lead  coordinates  in  that 
orientation.  This  is  a  simple,  but  conceptually  necessary,  step.  As  a  result,  the  model  conforms  to  the 
methods  used  by  M.  Coon  to  analyze  SAR  imagery.  He  has  shown  that  only  lead  coordinates  allow 
the  areas  to  be  calculated  accurately.  Of  course,  since  stretching  along  a  lead  must  be  zero  and  the 
trace  operator  is  invariant  under  coordinate  rotation,  the  rate  of  opening  or  closing  across  a  lead 
dj  must  equal  the  rate  of  dilating  tr  P; ,  expressed  in  either  lead  or  global  coordinates. 

Numerical  Method.  The  momentum  balance  equation  describes  the  sum  of  inertial,  Coriolis,  and  sea 
surface  tilt  accelerations  that  result  from  the  applied  air  and  water  stresses  and  divergence  of  the  ice 
stress  resultant. 

m(u  +  yk  x  U  -  gVr|)=  ofl  -  ow  +  V  6 

where  U  is  displacement,  U  is  velocity,  and  U  is  acceleration.  The  elastic  plastic  constitutive  law  is 

6  =  M^e 

where  Me/,  is  the  elastic  plastic  modulus  tensor.  This  modulus  tensor  results  from  combining  the  linear 
elastic  closure,  the  kinematic  relationship,  the  flow  rule,  and  yield  surface.  An  estimate  of  the  stress 
state  at  time  t+QAt  comes  from  a  difference  approximation  to  this  constitutive  law. 

The  Wilson  0-method  was  developed  as  an  implicit  finite  element  integration  scheme  for  dynamic 
structural  mechanics  problems.  It  may  be  derived  by  assuming  that  the  acceleration  is  linear  over  the 
time  period  ( t ,  t  +  0A/) 

(+TU-'U  =  —  (,+0A'U-fu) 

QAt v  ’ 

Integrating  the  linear  acceleration  provides  the  velocity  and  displacement  variable  as  a  function  of 
initial  values  and  the  acceleration.  Velocity  varies  as  a  quadratic  function  of  time 

f+Tu='u+'ux  +  — -  (,+0A'  tr-'u) 

20At v  ’ 


and  displacement  varies  as  a  cubic  function  of  time 


f+xu='u+'ux  +-fti  +  — -  (t+0At  u-fu) 

2  60At v  r 

The  velocity  at  time  t  +  0A t  is 

(+9ArU=f  U+'U0Af  +  —  ((+0ArU+ru) 

2  V  r 

and  the  displacement  is 

f+0ArU=  U+'U0At  +— rU  +  (f+9ArU-fu) 

2  6  v  / 

Finally,  we  solve  for  acceleration  and  velocity  in  terms  of  the  displacement  (all  at  time  t+QAt)  and 
known  values  at  time  t,  we  find  that 

'+0AfU  =  — (,+9A,U-'u)-— ' rU-2fU 
02Ar  v  ’  0At 

and 

f+eAfU  =  —  (f+eAfU-fu)-2'U  fU. 

0Af v  ’  2 

When  the  acceleration  and  velocity  at  time  f+0Af  are  expressed  in  terms  of  the  displacement,  the 
momentum  balance  equation  at  time  t+QAt  therefore  depends  only  on  the  displacement  (and  the  stress 
state) 

Momentum  balance  in  finite  element  method  notation  is 

MU  +  CU  +  KU  =  R 

where  M  is  the  generalized  mass  matrix,  C  is  the  generalized  damping  matrix,  and  K  is  the  generalized 
stiffness  matrix.  The  load  vector  on  the  right  hand  side  R  must  be  similarly  evaluated  at  this  time.  The 
displacement  at  time  t+QAt  can  be  determined.  In  practice,  we  can  reformulate  the  difference 
equations  in  terms  of  relative  displacement  over  the  time  step  0A t.  This  is  an  effective  approach  that 
was  used  by  Pritchard,  et  al.  [1983] 

If  the  parameter  0  exceeds  1.37,  the  scheme  is  unconditionally  stable.  This  property  is  important  for 
the  PIPS  3.0  ice  dynamics  model.  I  envisage  a  3-hour  time  step,  which  could  describe  all  essential 
processes  including  tidal  and  inertial  oscillations.  This  is  a  practicable  scheme  where  we  must  solve 
the  momentum  balance  equation  for  the  displacement  at  the  advanced  time.  The  method  is  therefore 
implicit  because  the  stress  state  must  be  calculated  at  the  same  time. 

IMPACT/APPLICATION 

Treating  sea  ice  as  an  anisotropic  material  provides  explicit  information  on  whether  or  not  leads  exist 
in  a  region,  and  their  orientations.  Knowing  when  and  where  leads  will  appear  is  useful  to  Navy 
operations.  Accounting  for  leads  (and  ridges)  explicitly,  will  allow  us  to  describe  ice  behavior  more 


accurately.  Perhaps  the  improved  accuracy  will  be  most  noticeable  in  lower  resolution  models  because 
the  discontinuous  behavior  will  be  better  described  within  each  grid  cell. 


TRANSITIONS 

The  constitutive  law  is  available  for  use  in  other  research  codes  and  in  Navy  operational  ice  forecasting 
systems.  I  will  work  with  Navy  personnel  and  contractors  to  convert  the  research  code  into  coding  that 
is  compatible  with  existing  codes  such  as  PIPS.  The  National  Ice  Center  has  identified  this  new 
anisotropic  constitutive  law  as  the  leading  candidate  for  including  in  PIPS  3.0. 

RELATED  PROJECTS 

The  technical  approach  requires  collaboration  with  M.  Coon  (Northwest  Research  Associates,  Inc.). 

A  formal  Management  Plan  was  prepared  to  describe  how  we  will  collaborate,  and  how  we  will 
prevent  duplication  of  our  efforts  so  that  we  meet  our  goals  with  the  scarce  resources  available. 

I  am  working  with  R.  Bourke,  J.  Wilson,  and  W.  Hibler  to  install  an  ambient  noise  model  into  future 
versions  of  PIPS.  I  am  helping  to  define  the  best  proxy  variables  available  in  the  PIPS  2  model. 
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